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Introduction 

Over  the  years,  computational  fluid  dynamics  (CFD)  has  steadily  gained  acceptance 
as  a  benefleial  tool  for  reducing  the  cost  of  designing  aerospace  vehicles.  As  the  accep¬ 
tance  of  CFD  grows,  many  new  applications  are  conceived  which  require  new  and  more 
sophisticated  algorithms  to  model  the  various  physical  processes  associated  with  advanced 
vehicle  design.  One  area  of  increasing  interest  is  the  design  of  hypersonic  vehicles,  such 
as  the  National  Aerospace  Plane  (NASP),  for  which  ground  based  wind  tunnel  testing  is 
either  technically  or  flscally  impractical. 

During  the  past  two  decades,  powerful  numerical  methods  for  solving  the  Navier- 
Stokes  equations  with  real  gas  effects  on  structured  meshes  have  been  developed  [1-2]. 
However,  during  the  same  time  period,  little  in  the  way  of  progress  has  been  accomplished 
in  the  area  of  efficiently  solving  the  Navier- Stokes  equations  with  real  gas  effects  on  unstruc¬ 
tured  meshes,  in  either  two  dimensions  or  three  dimensions.  Most  of  the  recent  successes  in 
unstructured  technology  have  been  in  the  area  of  solving  three  dimensional  inviscid  perfect 
gas  flows  or  two  dimensional  perfect  gas  viscous  flows.  While  technically  and  theoretically 
simple,  the  implementation  of  the  viscous  fluxes  along  with  finite  rate  chemistry  models 
in  practice  becomes  complicated  by  the  nature  of  unstructured  flow  solvers. 

For  the  work  presented  in  this  paper,  a  finite-volume  technique  has  been  chosen  in 
keeping  with  a  growing  trend  in  the  CFD  community.  Since  the  finite  volume  approach 
is  derived  from  the  integral  conservation  equations,  it  has  the  ability  to  properly  capture 
flow  discontinuities.  The  structure  of  this  particular  code  allows  the  computational  domain 
to  be  discretized  into  either  triangles,  quadrilaterals,  or  a  combination  of  both  in  two 
dimensions,  and  tetrahedra,  hexahedra,  prisms,  or  any  combination  of  these  elements  in 
three  dimensions. 

The  emphasis  of  the  current  work  centers  on  solution  algorithm  development.  While 
grid  generation  is  currently  one  of  the  limiting  factors  in  the  acceptance  of  unstructured 
flow  solvers  in  CFD,  suitable  techniques  are  becoming  more  readily  available. 

Mathematical  &;  Numerical  Formulation 

We  begin  the  development  of  our  numerical  algorithm  with  the  integral  form  of  the 
three  dimensional  Navier  Stokes  equations, 

where  Q  represents  the  vector  of  dependent  variables,  is  the  source  term  due  to  chemical 
reactions  and  non-equilibrium  thermodynamics,  F  •  ti  and  G  •  ti  represent  the  inviscid  and 
viscous  flux  of  mass,  momentum,  and  energy  out  of  the  control  volume  V  through  the 
surface  S  with  n  outward  unit  normal  from  S,  In  order  to  include  real  gas  effects,  the 
standard  system  of  five  equations  has  been  augmented  to  include  the  effects  of  a  generalized 
thermo-chemical  model  due  to  Grossman  et.  al.  [3].  For  our  purposes,  the  effect  is  to 
replace  the  global  continuity  equations  by  N  species  continuity  equations,  where  N  is  the 
number  of  species  in  the  chosen  chemistry  model.  In  addition,  to  account  for  M  species 
considered  to  be  in  thermal  non-equilibrium,  the  system  has  been  augmented  to  include 
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M  non-equilibrium  energy  conservation  equations.  In  light  of  this,  the  vector  of  dependent 
variables  and  thermo-chemical  source  terms  are  given  by; 


where  pi  is  the  species  density,  p  is  the  mixture  density,  u,  u,  and  w  are  the  components 
of  the  velocity  vector  in  the  principal  directions  z,  j,  and  k  respectively,  represents 
the  z*^  species  non-equilibrium  energy  contribution,  to  is  the  total  energy,  iVi  represents 
the  source  terms  due  to  chemical  reactions,  and  Qi  represents  the  source  terms  due  to 
vibrational  non-equilibrium  thermodynamics.  In  equation  (1)  we  assume  body  forces  are 
negligible.  The  governing  equations  are  closed  by  an  equation  of  state  relating  the  species 
densities,  pressure  and  energy^ 

The  inviscid  flux  vector  F  consists  of  components  from  the  three  principal  directions, 

F  =  fi  +  gj  -1-  hk  (3) 

where  the  flux  components  /,  g,  and  h  are  given  by: 


□  □ 


perfect  gas  flows  [4]  [5],  and  have  been  extended  to  flows  in  chemical  equilibrium  [6-9], 
and  to  flows  in  chemical  and  thermodynamic  non-equilibrium  [10-14].  Flux- vector  split 
techniques  have  proven  to  be  very  accurate  and  robust  when  used  for  hypersonic  flows  [15] 
and  are  fully  compatible  with  the  conservative  finite- volume  approach. 

Two  of  the  more  popular  flux- vector  split  algorithms  have  been  implemented  into  the 
unstructured  flow  solver.  These  are  the  algorithms  developed  by  Steger  and  Warming  [16] 
and  Van  Leer  [17].  The  basic  idea  of  both  methods  is  to  split  the  inviscid  flux  vector 
in  one  dimension  into  two  parts,  one  which  contains  the  information  that  propagates 
downstream  and  one  which  contains  the  information  that  propagates  upstream.  The  two 
parts  are  then  constructed  using  extrapolation  formulas  consistent  with  the  direction  of 
propagation.  Extension  of  the  algorithms  to  more  spatial  dimensions  is  accomplished  by 
superimposing  pseudo-one  dimensional  problems. 

In  addition  to  the  flux-vector  splitting  techniques,  the  flux  difference  splitting  tech¬ 
nique  originally  developed  by  Roe  [18]  has  been  implemented  into  the  unstructured  flow 
solver.  While  this  method  has  proven  to  be  less  robust  for  hypersonic  flows  [15],  it  has  been 
shown  to  be  more  accurate  than  the  flux-vector  split  techniques.  The  flux  difference  split¬ 
ting  technique  consist  of  an  approximate  Riemann  solver,  where  an  arbitrary  discontinuity 
is  supposed  to  exist  between  the  left  and  right  state.  An  approximate  solution  is  written 
for  this  situation  in  terms  of  waves  propagating  upstream  and  downstream.  In  this  case, 
the  flux  is  not  split,  but  reconstructed  from  the  upstream  and  downstream  contributions 
that  constitute  the  left  and  right  state  of  the  Riemann  problem. 

Similar  to  the  inviscid  flux  vector,  the  viscous  flux  vector  G  may  be  written  in  com¬ 
ponent  form  as: 

<?  =  /«*  +  9vj  +  Kk  (5) 

where  the  viscous  flux  components  f„,  gv,  and  hy  are  given  by: 


In  this  formulation,  mass  diffusion  has  been  neglected.  The  stresses  in  equation  (6)  are 
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written  in  terms  of  the  velocity  gradients  via: 

Txx  =  A*  (^Ux  —  -  {Ux  +  Vy  + 

r„  =  ^  (2V,  -l(U.  +  Vy  +  W.)^  (6a) 

T2Z  =  —  -  (Ux  +  Vy  + 

Txy  =  f^{Uy  +  Vx) 

T,,  =  fiiU,  +  W,)  (66) 

Ty,  =(,{¥.  +  Wy) 

Tx  —  UTxx  “1"  1^ "^xj/  d"  Id^'^xx 

fy  =  UTxy  -\-VTyy+  WTj,^  (6c) 

=  C/Tiz  +  VTy:  +  WTzz 

The  heat  flux  terms  in  equation  (6)  are  written  in  terms  of  the  temperature  gradients  via: 

Qx  =  -kTx 

q,  =  -kT,  (6<J) 

qz  =  -kTz 

where  fi  is  the  laminar  viscosity,  k  is  the  thermal  conductivity,  and  T  is  the  translational 
temperature.  The  subscript  x,  y,  and  z  represent  the  gradient  of  the  term  in  the  corre¬ 
sponding  spatial  direction.  The  bulk  viscosity,  A  =  — 2^/3  is  based  on  Stokes  hypothesis, 
which  is  valid  when  Newtonian  fluids  are  considered  and  bulk  viscosity  effects  are  neglected. 
Bulk  viscosity  effects  should  be  accounted  for  when  rotational  non-equilibrium  is  present 
[19].  In  addition  bulk  viscosity  effects  exhibit  a  significant  influence  on  sound  propagation 
and  shock- wave  structure  [19]. 

Spatial  Discretization 

Since  the  finite- volume  method  is  simply  a  discrete  representation  of  the  integral  con¬ 
servation  laws,  it  has  the  property  of  properly  capturing  all  flow  discontinuities.  In  the 
finite- volume  formulation,  the  domain  is  subdivided  into  a  discrete  number  of  control  vol¬ 
umes.  The  variation  of  the  flow  quantities  may  then  be  approximated  by  a  polynomial  over 
each  control  volume.  In  this  work  we  use  a  polynomial  that  satisfies  several  reconstruction 
criteria  including  k- exactness.  The  discrete  equations  are  then  solved  simultaneously  over 
each  individual  control  volume.  The  particular  functions  describing  the  flow  quantities 
within  each  control  volume  are  used  to  evaluate  the  volume  and  surface  integrals  of  the 
governing  equations.  The  choice  of  these  functions  determine  the  spatial  accuracy  of  the 
solution. 

In  discretizing  the  governing  equations,  the  conserved  variables  Q  must  be  integrated 
over  each  control  volume,  the  fluxes  must  be  integrated  over  all  of  the  faces  surrounding 
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each  control  volume,  and  the  source  terms  W  must  be  integrated  over  each  control  volume. 
The  inviscid  fluxes  are  discretized  using  a  characteristic  based  flux-split  algorithm,  either 
Steger- Warming,  Van  Leer,  or  Roe’s  numerical  flux  function.  All  of  these  algorithms 
require  that  the  dependent  variables  be  interpolated  from  either  side  of  each  cell  face. 
Therefore,  procedures  need  to  be  developed  which  model  the  point-wise  spatial  variation 
of  conservative  variables  Q  in  terms  of  the  known  cell  averaged  conservative  variables  Q. 
The  key  to  the  success  of  any  finite  volume  flow  solver  is  the  ability  to  accurately  produce 
these  point-wise  variations  in  the  flow  variables. 


Reconstruction 

The  problem  of  determining  the  point-wise  variation  of  the  flow  variables  from  the 
cell  averages  is  known  as  the  reconstruction  problem.  Simply  stated,  the  reconstruction 
problem  involves  the  process  of  determining  the  point-wise  distribution  given  the  cell  aver¬ 
ages  with  the  restriction  that  integrating  the  point- wise  function  recovers  the  cell  averages. 
Mathematically  this  restriction  becomes  one  of  ensuring  conservation  of  the  mean 

Q  =  ^  JJJ^Qix,y,z)  dV  (7) 

For  steady  state  solutions,  the  accuracy  of  the  scheme  is  determined  by  the  accuracy  of 
the  expressions  used  to  evaluate  the  fluxes  at  the  cell  interface. 

In  the  unstructured  flow  solver,  three  higher  order  reconstruction  jnethods  are  em¬ 
ployed.  The  first  two  are  gradient  based  schemes  proposed  by  Barth[20]  and  Frink[21].  For 
the  problems  presented  here,  those  methods  were  passed  up  in  favor  of  a  k-ex&ci  formula¬ 
tion  originally  proposed  by  Barth[22].  The  flavor  of  the  fc-exact  formulation  implemented 
in  the  unstructured  code  was  developed  by  Mitchell  and  Walters  [23].  For  this  method, 
the  functional  form  for  Q{x,y,z)  is  a  polynomial  expression  of  degree  1  for  second  order 
accuracy. 

Q{x,y,z)  =Co  +  Cix  +C2y  +  C3Z  (8) 

The  use  of  the  polynomial  expression  in  equation  (8)  has  three  advantages  which  make  it 
a  desirable  choice  for  implementation.  First,  the  expressions  are  easily  integrated  over  a 
control  volume  via  simple  quadrature  formulas.  Second,  the  computation  of  weak  solutions 
is  permissible  via  two  independent  states,  Ql  and  Qr,  which  are  evaluated  at  the  cell 
interface.  Discontinuities  may  then  be  resolved  using  an  appropriate  Riemann  solver.  Last 
and  perhaps  most  important  for  our  discussion,  they  are  easily  differentiated  with  respect 
to  the  spatial  coordinates  x,  y,  and  2.  This  allows  the  computation  of  the  gradients  in  the 
viscous  fluxes  to  be  performed  with  little  additional  computational  effort. 

In  equation  (8),  there  are  four  coefficients  for  each  polynomial  which  must  be  solved 
for  (in  two  dimensions  the  number  of  coefficients  is  reduced  to  three).  Therefore,  we  need 
four  cells  in  which  we  can  apply  the  restriction,  equation  (7),  to  arrive  at  a  linear  system 
for  the  coefficients.  For  the  second  order  accurate,  k  =  I  reconstruction  considered  here. 


5 


the  four  restrictions  on  our  reconstruction  polynomial  becomes; 


1 

Vi 

1 

V2 

1 

1 


x,y,z)  dx  dy  dz  =  Qi 
x,y,z)  dx  dy  dz  =  Q2 
x,y,z)  dx  dy  dz  =  Q3 
x,y,z)  dx  dy  dz  = 


(9) 


where  the  subscript  1  denotes  the  first  cell  in  the  stencil,  2  the  second  cell  in  the  stencil, 
and  so  on.  In  equation  (9),  the  polynomial  function  for  Q{x,  y,  z)  is  the  same  for  each  of  the 
equations.  For  our  purposes,  any  acceptable  second  order  quadrature  rule  is  satisfactory. 


^  JJJ^Qix,y,z)  =  Q{x,y,z)  +  0{h'^)  (10) 

where  x,  y,  and  z  specify  the  location  of  the  cell  centers.  Therefore,  when  the  quadra¬ 
ture  rule  is  applied  to  equation  (9)  we  arrive  at  a  linear  system  for  the  reconstruction 
polynomials. 

r  1  xi  yi  zi 
1  X2  m  Z2 
1  X3  yz  zz 
.1  Xi  yi  Zi 

Equation  (11)  can  be  solved  directly  via  symbolic  inversion.  Substitution  into  equation  (8) 
yields  an  expression  for  Q{x,  y,  z)  in  terms  of  the  centroids  of  the  stencil  cells  and  the  cell 
average  variables.  For  a  stencil  which  is  fixed  in  time,  the  only  part  of  the  expression  which 
changes  are  the  values  of  the  cell  averaged  variables.  Therefore,  we  may  rewrite  equation 
(8)  as  a  vector  of  reconstruction  weights  u;  and  a  vector  of  the  cell  averaged  variables  Q: 


(11) 


u?i(x,y,z)/ 

where  the  spatial  dependency  come  in  the  form  of  the  reconstruction  weights.  With  a 
simple  second  order  accurate  flux  integration,  the  weights  may  be  determined  once  at  the 
beginning  of  a  solution  thereby  reducing  the  computational  effort  and  storage  requirements. 


Discretization  of  the  Viscous  Fluxes 

The  discretization  of  the  viscous  fluxes  is  dependent  upon  the  method  of  reconstruc¬ 
tion  used  to  determine  the  point-wise  variation  of  the  conservative  variables  from  the 
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cell  averaged  values.  The  discussion  in  this  section  is  based  on  the  k-exact  reconstruc¬ 
tion  method,  however,  the  methods  developed  are  general  enough  that  extension  to  other 
reconstruction  methods  is  straightforward. 

For  the  k-exact  reconstruction  method,  the  gradients  may  be  determined  by  simple 
differentiation  of  the  reconstruction  polynomials  with  respect  to  the  spatial  coordinates. 
In  terms  of  the  weights,  the  gradients  of  the  u  velocity  component  become: 


Ux  =  +(^2^U2  +I^3^U3  +U)4^U4 

Uy  =  +  U!2yU2  +  (jJ3yU3 (13) 

Uz  —  1^1  z  Ui  i*}2z^2  +  ^3z^3  + 

where  is  the  partial  derivative  of  the  reconstruction  weight  with  respect  to  the  x  spatial 
coordinate.  The  other  gradients  in  equations  6(a)-6(d)  are  determined  in  a  similar  fashion. 

The  transport  properties  are  computed  from  the  reconstructed  values  at  the  interface. 
In  other  words,  for  Sutherland’s  law,  the  temperature  at  the  cell  interace  is  first  computed, 
from  which  the  viscosity  and  thermal  conductivity  may  then  be  computed.  This  eliminates 
the  need  to  reconstruct  the  transport  properties  themselves.  The  computation  of  the 
species  and  mixture  transport  properties  are  discussed  shortly. 

Once  the  gradients  are  computed  the  viscous  fluxes  may  be  computed  via: 


/ 


0 

0 

0 


TjxTxx  "b  f}yTxy  "t"  fjz'^xz 
f)xTxy  +  Vy'^yy  "b  Vz'^yz 
^x'^xz  *b  'Hy'^yz  ”b  '^z'^zz 
0 
0 

0 

-  9.)  +  %  (Tj  -  ly)  +  <1y  (t,  -  qy)  / 


(14) 


where  rj  is  the  outward  normal  and  rjx  is  the  x-component  of  the  outward  normal. 


Viscosity  Coefficient 

The  laminar  viscosity  coefficient  (/i)  is  computed  for  each  species  in  the  chemistry 
model  using  curve  fits  which  are  based  on  experimental  data.  The  advantage  of  using  curve 
fits  is  their  simplicity  and  availability  of  semi-empirical  rules  for  recovering  the  mixture 
values.  Two  experimental  curve  fits  are  available  in  the  unstructured  code.  The  first  is 
one  of  the  most  widely  adopted  curve  fits  due  to  Blottner  [24], 

fx,  z=  exp[{As\nT  +  Bs)\nT  +  Cs]  s  =  l,...,JV  (15) 
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where  Ag,  Bs,  and  Cs  are  coefficients  determined  from  fitting  the  experimental  data. 
Another  popular  experimental  curve  fit  is  due  to  Sutherland, 


Ms  =  T 


1.5 


Es 


T  +  F/ 


(16) 


where  Eg  and  Eg  are  empirically  derived. 

The  mixture  viscosity  coefficient  is  determined  via  Wilke’s  semi-empirical  rule  [25], 
which  is  based  on  kinetic  theory, 


(17) 


where 


SJ 


and  Xg  is  the  mole  fraction  of  species  s. 


(18) 


Thermal  Conductivity  Coefficient 

Like  the  species  viscosity,  the  species  thermal  conductivity  is  computed  from  experi¬ 
mental  curve  fits.  One  such  curve  fit  is  Eucken’s  relation. 


k.  =  ,  ‘>  =  h',N  (19) 

where  ^  is  the  translational  contribution  to  the  specific  heat  at  constant  volume.  For  a 
perfect  gas,  with  7  =  7/5,  Eucken’s  formula  returns  a  Prandtl  number  Pr  =  0.737,  which 
is  in  good  agreement  with  experimental  values. 

A  second  curve  fit,  due  to  Sutherland,  has  been  implemented  for  determining  the 
species  thermal  conductivity.  The  form  of  this  curve  fit  is  given  by 


h  - 

T  +  Fg' 


s  =  l,...,iV 


(20) 


where  the  constants  Eg  and  Eg  are  determined  from  experimental  data. 

The  mixture  thermal  conductivity  is  computed  using  Wilke’s  mixture  rule  [25], 


N 


S=1 


(21) 


where  the  (f^g^j  are  the  same  as  in  the  formulation  for  the  mixture  viscosity  coefficient. 

A  simplified  method  for  determining  the  mixture  thermal  conductivity  has  been  im¬ 
plemented  by  assuming  a  constant  Prandtl  number  and  using  the  definition  of  the  Prandtl 
number, 

k  =  where  Pr  =  constant  (22) 

Pr 
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In  this  formulation,  the  species  thermal  conductivities  are  not  required. 


Stencil  selection 

One  of  the  more  challenging  tasks  associated  with  the  k-exact  formulation  is  the 
selection  of  the  cells  which  complete  each  stencil.  For  the  inviscid  fluxes,  two  such  stencils 
are  required,  one  for  the  left  state  and  one  for  the  right  state.  For  the  viscous  stencil,  only 
one  stencil  is  required,  typically  a  centered  stencil. 

Figure  1(a)  shows  a  two  dimensional  triangular  mesh  with  a  representative  upwind 
stencil.  In  flgure  1(a),  the  bold  face  represents  the  face  for  which  the  stencil  must  be  gen¬ 
erated.  The  circled  numbers  represent  the  cells  which  would  complete  the  two  dimensional 
stencil.  In  two  dimensions  three  cells  are  required  to  complete  the  system  of  equations  for  a 
fc  =  1  reconstruction.  In  three  dimensions  four  cells  are  required  for  a  =  1  reconstruction. 

The  inviscid  stencil  selection  process  begins  by  gathering  the  cells  adjacent  to  the  face, 
in  this  case,  cell  Ci  and  C2  which  we  refer  to  as  the  parent  cell  and  biased  cell  respectively. 
The  next  step  is  to  gather  the  cells  surrounding  the  parent  cell  {Cu  and  Cu)-  For  the 
case  of  a  triangular  mesh,  the  selection  process  would  simply  choose  the  parent  cell  and 
its  two  surrounding  cells.  The  stencil  would  then  be  tested  to  ensure  it  was  a  valid  stencil 
and  if  it  failed  either  the  biased  cell  would  be  used  or  the  surrounding  cells  of  Cn  and 
C12  would  be  added  to  the  possible  list  of  cells.  For  quadrilateral  meshes,  where  there  are 
four  cells  from  which  to  choose,  the  dot  product  of  the  cell  centroid  with  the  face  metric  is 
used  as  a  basis  for  determining  which  cells  to  include  in  the  stencil.  In  this  case,  the  two 
surrounding  cells  with  the  smallest  dot  product  and  the  parent  cell  is  typically  used.  In 
three  dimensions  the  process  becomes  complicated  by  the  fact  that  four  cells  are  required 
and  poor  stencils  show  up  more  readily  than  in  two  dimensions.  In  the  unstructured  code, 
several  different  methods  for  determining  the  inviscid  and  viscous  stencil  are  provided.  In 
three  dimensions,  it  was  found  that  a  condition  number  test  on  the  4x4  coefficient  matrix 
produced  the  best  results  to  date.  This  requires  an  extra  step  in  the  stencil  selection 
process.  After  gathering  all  of  the  cells  which  are  to  be  considered  for  inclusion  into  the 
stencil,  we  need  to  construct  a  group  of  valid  stencils  for  which  we  are  going  to  apply  the 
condition  number  test.  Instead  of  testing  every  possible  stencil  which  can  be  an  enormous 
task,  the  algorithm  can  be  simplified  by  restricting  how  we  define  our  stencils.  In  our  case, 
we  consider  only  stencils  made  up  of  cells  which  share  a  face  with  at  least  one  other  cell  in 
the  stencil.  We  use  the  term  continuous  stencil  to  refer  to  these  types  of  stencils.  At  this 
point  it  is  too  early  to  say  whether  a  continuous  stencil  provides  a  more  stable  solution 
algorithm. 

Figure  1(b)  shows  a  two  dimensional  triangular  mesh  with  a  representative  viscous 
stencil.  In  this  figure,  the  bold  face  represents  the  face  for  which  the  flux  is  being  computed. 
The  circled  numbers  represent  the  cells  which  would  be  chosen  to  complete  the  stencil. 

For  the  viscous  stencil,  a  centered  stencil  is  desirable.  The  selection  process  is  similar 
to  the  inviscid  stencil  with  the  exception  that  the  cells  surrounding  both  the  parent  and 
biased  cell  are  gathered.  The  stencil  is  then  chosen  from  the  parent,  biased,  and  one  of 
the  surrounding  cells.  As  with  the  inviscid  stencil  selection,  test  may  be  applied  directly 
to  the  cells  or  to  the  stencils. 
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Chemical  Source  Terms 

Following  the  work  of  Grossman  and  Cinnella  [3],  we  have  implemented  thermo¬ 
chemical  models  into  the  flow  solver.  A  brief  description  of  the  models  is  presented  here 
for  completeness.  A  fairly  general  simulation  of  chemical  effects  can  be  included  in  the 
governing  equations,  for  a  system  containing  N  species  in  which  J  reactions  tahe  place 

H - 1-  ^N,j^N  ^ 

I'ljXi  -t-  V2,j^2  H - H  ,  (23) 

j  —  I5  2,  •  •  •  ?  J) 


The  i/-  •  and  the  u" j  are  the  stoichiometric  coefficients  of  the  reactants  and  products  of 
species  JAj  in  the  jth  reaction,  respectively.  In  the  source  vector  W ,  the  rate  of  production 
of  species  i  may  be  written  as  [26] 


i=i 


''i.j 


i  =  1,2,. ..,iV 


(24) 


where  for  reaction  j  the  forward  and  backward  reaction  rates,  kfj  and  are  assumed 
to  be  known  functions  of  temperature,  and  are  related  by  thermodynamics 


kb,j 


(25) 


The  equilibrium  constant, is  a  known  function  of  the  thermodynamic  state.  Using 
this  relation,  we  obtain 


i=i 


''i.j 


i  =  1,2,. ..,iV 


(26) 


Note,  the  term  in  the  brackets  goes  to  zero  when  equilibrium  compositions  are  reached  for 
Pi,  since  it  reduces  to  a  formulation  of  the  Law  of  Mass  Action 


KeJ  = 


nil, 

n£,  ’ 


(27) 
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valid  for  equilibrium  concentrations. 

The  unstructured  flow  solver  currently  contains  28  chemistry  models,  of  which  5  are 
air  chemistry  and  18  are  hydrogen-air  chemistry.  In  addition,  there  are  approximately  450 
reactions  and  34  species  from  which  to  choose. 

Internal  Energy 

Assuming  each  species  in  the  gas  mixture  behaves  as  a  thermally  perfect  gas,  the 
internal  energy,  ti  may  be  defined  as 


e,-  =  e.(T)  4-  e„,  (28) 

where  is  the  known  equilibrium  portion  of  the  energy,  T  is  the  translational  temperature, 
and  is  the  non-equilibrium  portion  of  the  energy.  The  portion  of  the  internal  energy 
due  to  electronic  excitation  has  been  neglected  [19]. 

It  is  convenient  to  express  the  equilibrium  portion  of  the  internal  energy,  ei ,  in  terms 
of  a  specific  heat  at  constant  volume,  c„,. ,  as  follows 


U  =  f  c„.(r) 


dr  +  /t/i 


(29) 


where  —  dci/dT^  and  hf^  is  the  heat  of  formation  of  species  i, 

The'internal  energy  of  a  gas  mixture  composed  of  N  species,  with  the  first  M  species 
assumed  to  contain  a  non-equilibrium  portion  of  their  internal  energy  may  be  written  as 


e 


N 


N 


M 


=  E-'<+Et'- 

i=l  ^  2=1  ^ 


The  mixture  gas  constant  is  given  by 


R 


(30) 


(31) 


Equation  of  State 

A  very  important  assumption  that  has  been  made  is  the  consideration  of  weakly  in¬ 
teracting  particles,  where  the  interaction  potential  is  decaying  very  rapidly  in  space  and  is 
negligible  after  a  typical  distance  of  the  order  of  the  particle  radius.  The  most  noticeable 
effect  of  this  assumption  is  that  the  mixture  components  will  behave  as  thermally  perfect 
gases,  where  the  internal  energy  is  only  a  function  of  temperature.  Strictly  related  to 
this  result  is  the  Aralidity  of  Dalton’s  Law ,  whereby  the  mixture  pressure  is  the  summation 
of  partial  pressures.  The  applicability  of  the  weak-interaction  assumption  is  restricted  to 
conditions  of  low  density  and  moderate  to  high  temperature.  Many  flows  of  engineering 
interest  meet  these  requirements. 
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(32) 


Consequently,  the  relationship  between  pressure  and  temperature  will  read 

N 

p  =  ^  PiRiT  =  pRT 

i=l 


where  the  mixture  density  is  given  by 


P 


(33) 


and  where  the  mixture  gas  constant,  R,  is  defined  by  (31).  In  equation  (32),  we  have  made 
the  simplifying  assumption  that  the  translational  temperature  is  the  same  for  each  species, 
Ti  =  T.  This  simplification  is  reasonable  for  flows  where  the  mass  of  the  molecules  in  the 
mixture  are  on  the  same  order  of  magnitude.  Thus,  this  will  not  hold  for  flows  with  free 
electrons. 

The  state  relationship  of  the  pressure  to  the  specific  internal  energy  occurs  implicitly 
through  the  temperature.  For  a  given  chemical  composition,  internal  energy  and  non¬ 
equilibrium  energy,  the  temperature  must  be  evaluated  from 


N 


M 


(34) 


ST  f 


i=l 


Iterative  procedures  could  then  be  used  to  solve  for  T .  Once  T  is  found,  the  pressure  may 
be  evaluated  from  (32).  However,  our  unstructured  code  integrates  the  primitive  variables 
Pi,  u,  u,  w,  e„( ,  and  p  in  time.  Therefore,  this  additional  iterative  procedure  is  not  required 
since  the  temperature  T  can  be  found  directly  from  the  pressure  p.  The  implementation  of 
time  integration  schemes  applied  to  the  primitive  variables  while  retaining  the  conservative 
property  is  straightforward. 


Non-Equilibrium  Thermodynamics 

The  nonequilibrium  thermodynamic  model  groups  the  translational  and  rotational 
contributions  in  ,  thus  considering  them  in  equilibrium  at  the  translational  temperature 

T. 

ti  =  UiRiT  +  hf.,  i  —  (35) 

In  this  equation,  the  constants  rii  are  given  by 


ni=  < 


for  monotomic  species, 

for  diatomic  &  linear  polyatomic  species, 

for  nonlinear  polyatomic  species. 


Vibrational  contributions  are  assumed  to  be  in  non-equilibrium  and  electron  excitation 
has  been  neglected.  As  a  result,  this  model  should  be  restricted  to  use  with  air  chemistry 
models. 
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The  elastic  contribution  to  the  source  terms  in  the  non-equilibrium  vibrational  en¬ 
ergy  equation  can  be  modeled  by  considering  the  exchanges  between  the  vibrational  and 
translational-rotational  contributions  to  the  internal  energy.  The  rate  of  energy  exchange 
between  the  vibrational  and  translational-rotational  modes  are  assumed  to  be  of  a  Landau- 
Teller  form, 

,  s  =  (36) 

/  n,E  T~s 

where  e*  is  the  equilibrium  vibrational  contribution  to  the  internal  energy  at  the  transla¬ 
tional  temperature  T 

e;  =  (37) 

and  r  is  a  relaxation  time  given  by  Millikan  and  White  [27] 


r.  = 


i.omo® 

p 


exp 


-  18.42 


(38) 


where 


=  1.16  Pss  = 


MsMs 

Ms  + 


(39) 


The  inelastic  contribution  can  be  modeled  by  taking  into  account  chemical  reactions 


(^*)n 

where  no  attempt  has  been  made  to  consider  coupling  effects  though  equivalent  tempera¬ 
tures  or  by  other  means.  The  total  source  term  for  the  non-equilibrium  energy  equation 
is  given  by  the  sum  of  the  elastic  and  inelastic  contributions, 


(41) 


Equilibrium  Thermodynamics 

Two  thermal  equilibrium  models  are  implemented  in  the  unstructured  code.  The  first 
model  is  a  thermal  equilibrium  model  where  Cn,  =  0  and  vibrational  contributions  are 
included  into  e,-  by  means  of  a  simple  harmonic  oscillator  formula 


Ct  =  TiiRiT  -|- 


RiQv,i 


^  +  hf^ ,  i  —  1, . . .  5  iV 


(42) 


where  0t,,,  are  characteristic  vibrational  temperatures.  More  than  one  characteristic  tem¬ 
perature  will  appear  in  the  formulas  for  polyatomic  molecules.  In  this  model  we  have 
neglected  electronic  excitation.  Therefore,  this  model  is  applicable  to  air  chemistry  mod¬ 
els. 
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The  second  model  is  a  thermal  equilibrium  model  where  polynomials  are  used  to  curve 
fit  Cvi  [28] 

Cvi  =  0,1  +  02  IT  +  03  +  a^T  + 

+  oeT*  +  +  ogT^ 

from  which  the  internal  energy  can  be  found  from 


For  the  34  species  currently  implemented  in  our  database,  the  valid  temperature  range  for 
the  curve  fits  is  from  200  K  to  6,000  K. 

Turbulence  Modeling 

With  the  growing  acceptance  and  use  of  computational  fluid  dynamics,  many  new 
applications  are  conceived  which  require  new  and  more  sophisticated  algorithms  to  model 
the  various  physical  processes  associated  with  current  vehicle  design.  Since  most  of  the 
flows  of  interest  in  practice  are  turbulent,  one  component  required  to  successfully  model 
these  complex  flows  is  the  modeling  of  turbulence.  In  this  work,  we  have  implemented  the 
Spalart  and  Allmaras  one  equation  turbulence  model  [29]. 

One  equation  models  use  one  additional  field  equation  to  specify  a  velocity  scale. 
Historically  the  velocity  scale  chosen  for  one-equation  models  is  the  turbulent  kinetic  energy 
(TKE).  However,  one  equation  models  based  on  the  turbulent  kinetic  energy  still  require 
an  algebraic  relation  for  the  length  scale  which  leads  to  results  not  very  different  from  the 
algebraic  models.  Like  the  algebraic  models,  these  models  are  cumbersome  to  implement 
in  an  unstructured  format.  A  new  group  of  one  equation  models  has  been  developed  in 
which  the  field  equation  represents  a  transport  equation  for  the  turbulence  viscosity.  These 
new  one  equation  models  are  ^^locaP^  in  that  the  equation  at  one  point  does  not  depend 
on  the  solution  at  other  points  and  therefore  are  compatible  with  both  structured  and 
imstructured  grid  topologies. 

Transport  Equation 

A  summary  of  the  original  model  as  presented  in  the  original  Spalart  and  Allmaras 
paper  [29]  is  briefly  reviewed  here.  The  differential  transport  equation  along  with  all 
relevant  closure  coefficients  and  auxiliary  relations  are  provided.  Immediately  following 
the  discussion  on  the  differential  form  of  the  governing  equation  is  the  integral  form  of  the 
transport  equation  suitable  for  a  finite  volume  fluid  dynamics  code. 


Differential  Formulation 

Our  discussion  closely  follows  the  appendix  of  the  Spalart  and  Allmaras  paper  with  two 
exceptions.  The  implementation  followed  here  does  not  include  the  optional  trip  function 
required  for  transition.  Due  to  the  removal  of  this  term,  the  boundary  and  initial  conditions 
must  also  be  slightly  modified  from  those  discussed  in  the  appendix  of  the  original  paper. 
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The  one-equation  model  provides  a  single  differential  transport  equation  for  the  tur¬ 
bulent  quantity,  v.  The  dependent  variable  in  the  transport  equation  is  related  to  the  eddy 
viscosity  by:  ^ 

=  i>fvu  fvi  =  3  ^ 

X  I  ^vl 

where  i>  is  the  dependent  (working)  variable  and  u  is  the  molecular  viscosity  (fi/p).  The 
differential  form  of  the  transport  equation  for  the  dependent  variable  is  given  by: 


Du 

Dt 


=  CbiSu  -  C^ifu 


d 


+ 


1  r 


V  ((i^  +  i>)  V  ^)  +  Cb2 


(46) 


The  difference  between  equation  (46)  and  equation  (A2)  in  [29]  is  the  absence  of  the 
tripping  functions  required  for  modeling  transition.  This  is  equivalent  to  setting  the  closure 
coefficients  Cn  and  Cts  to  zero  in  the  equations  provided  in  the  appendix  of  [29] 

The  transport  equation  (46)  represents  on  the  left  side  a  local  time  rate  of  change 
and  convective  acceleration  via  the  material  derivative,  and  on  the  right  hand  side  from 
left  to  right  a  production  term,  destruction  term,  and  two  diffusion  terms,  respectively. 
The  auxiliary  relations  required  to  complete  the  production  term  are  given  by: 


5  =  5  + 


nu 


fv2 


fv2  =  1  ~ 


X 

1  +  Xfvl 


(47) 


where  5  is  the  magnitude  of  the  vorticity  and  given  in  terms  of  the  velocity  gradients  using 
tensor  notation  /  ott  arr  \ 

5  =  ^2nijnij,  = 

and  d  is  the  distance  to  the  nearest  solid  boundary. 

The  auxiliary  relations  required  to  complete  the  destruction  term  are  provided  by: 


fw  =  9 


1  +  Cu,3® 

J 


.1/6 


g  =  r  +  Cu,2  (r®  -  r)  , 


nu 


(49) 


For  large  values  of  r,  the  function  fw  reaches  a  constant.  Therefore,  values  of  r  greater 
than  10  may  be  truncated. 

The  model  is  completed  by  the  closure  coefficients.  The  coefficients  are  C(,l  =  0.1355, 
Cb2  =  0.622,  K  =  0.41,  <7  =  2/3,  c,„l  =  C61/k^  +  (1  +  C62)/<t,  Cy,2  =  0.3,  =  2,  and 

c„i  =  7.1. 

Turbulent  heat  transfer  is  computed  using  a  constant  turbulent  Prandtl  number  input 
by  the  user  (typically  equal  to  0.9).  The  turbulent  thermal  conductivity  is  then  computed 
using  the  definition  of  the  Prandtl  number 


kt  =  ^ 


(50) 
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where  kt  is  the  turbulent  thermal  conductivity,  fXt  =  pvt  is  the  turbulent  viscosity,  Cp  is 
the  mixture  specific  heat,  and  at  is  the  turbulent  Prandtl  number. 

The  boundary  conditions  for  the  turbulence  viscosity  are  simple  and  easily  imple¬ 
mented.  For  wall  bounded  flows,  the  wall  boundaxy  condition  is  to  set  the  turbulence 
viscosity  to  zero  {v  =  0).  All  other  boundary  conditions  axe  fixed  at  the  freestream.  The 
initial  condition  is  set  to  the  freestream  turbulence  viscosity.  The  freestream  turbulence 
viscosity  is  input  by  the  user,  however,  a  value  roughly  equal  to  10  times  the  molecular 
viscosity  has  been  found  to  produce  adequate  results.  It  should  be  noted  that  this  is  dif¬ 
ferent  from  the  values  proposed  in  the  original  model.  The  values  provided  in  the  original 
model  assume  the  triping  term  is  implemented. 

The  integral  conservation  form  of  the  turbulence  model  transport  equation  may  be 
written  as: 


Cwifu 


s  \<^ 


The  auxiliary  functions  and  closure  coefficients  remain  the  same  for  the  integral  transport 
equation. 


Convective  Flux 

The  convective  flux. 


has  been  implemented  into  the  unstructured  flow  solver  with  two  different  methods.  In 
both  methods,  regions  where  the  Mach  number  is  greater  than  unity,  a  full  upwind  biased 
flux  is  used.  In  regions  where  the  Mach  number  is  less  than  unity,  either  a  full  flux  is  used 
or  a  flux  vector  splitting  is  used  depending  on  user  preference.  At  this  point  it  is  too  early 
to  determine  whether  one  method  has  a  distinct  advantage  over  the  other. 

The  left  and  right  state  vectors  (solution  and  thermodynamic  property  arrays)  are 
determined  using  a  first  order  reconstruction.  If  we  denote  the  x,y,  and  2  metrics  by  rjx, 
rjy,  and  r]z  respectively,  we  compute  the  left  and  right  Mach  numbers  via: 

Ui  =  r)xUi  +  rjyVi  +  rizWi  Mi  =  Ui/ai 
Ur  =  'nxUr  +  'nyVr  +  rizWr  Mr  =  Urfo,r 

where  w,  v,  and  w  axe  the  x,  y,  and  2  components  of  velocity,  M  is  the  Mach  number,  and 
a  the  speed  of  sound.  Based  on  the  magnitude  and  sign  of  the  Mach  numbers,  the  method 
of  computing  the  convective  flux  is  determined. 

If  the  left  state  Mach  number  is  greater  than  unity  (M/  >  1)  then  the  flux  is  computed 
using  a  full  upwind  biased  stencil, 

F  =  Uivi  (54) 
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and  in  a  similar  fashion,  if  the  right  state  Mach  number  is  less  than  -1  then  the  flux  is 
computed  using  , 

F  =  UrVr  (55) 

Obviously,  in  both  cases  the  flow  is  supersonic  and  an  upwind  stencil  is  desirable. 

Regions  where  the  flow  is  subsonic,  the  convective  flux  is  computed  using  either  a  full 

flux, 

F  =  i  {Urn  +  UrVr)  (56) 

or  a  split  flux  which  follows  the  Van  Leer  flux  vector  splitting  scheme  discussed  in  the 
earlier  sections, 

F  =  Q  {Ml  +  1)"  a)j  VI  -  J  {Mr  -  1)"  ar)  Vr  (57) 


Diffusive  Flux 

The  diffusive  flux, 

JJ  ^ 

is  computed  in  a  similar  fashion  to  the  viscous  fluxes.  The  same  reconstruction  weights 
that  were  used  to  compute  the  viscous  fluxes  are  used  to  compute  the  turbulence  model 
diffusive  flux.  Recall,  that  in  general  there  are  four  reconstruction  weights  for  =  1 
reconstruction  in  three  dimensions.  The  variables  are  reconstructed  to  the  face  via  the 
reconstruction  weights, 

Vfact  =  iJJlVl  +  <jJ2V2  ‘^SVZ  +  <^iVi  (59) 

where  a;,-  is  the  reconstruction  weight  for  the  ith  cell  in  the  reconstruction  stencil  and 
Vi  is  the  cell  centered  value  of  the  reconstructed  variable.  The  molecular  viscosity  is 
reconstructed  in  the  same  manner. 

The  gradients  are  computed  using  the  partial  derivatives  of  the  reconstruction  weights. 
Recall  that  for  A:  =  1  and  greater  reconstructions,  the  weights  are  a  function  of  spatial 
coordinates  x,  y,  and  z.  Therefore,  the  gradients  of  the  turbulence  viscosity  may  be 
computed  using 

+  <-^2xii'2  +  ^31^3  +  (60) 

ox 

where  ujix  is  the  partial  derivative  of  the  ith  reconstruction  weight  with  respect  to  the  x 
spatial  coordinate.  The  gradients  in  the  other  spatial  coordinate  directions  are  obtained 
in  the  same  manner.  The  formulation  of  the  reconstruction  weights  and  their  derivatives 
for  a  A:  =  1  reconstruction  are  given  in  [30]. 


Distance  Function 

From  an  unstructured  coding  point  of  view,  the  computation  of  the  distance  function 
represents  the  greatest  challenge  in  implementing  the  one  equation  turbulence  model.  The 
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distance  function  may  be  computed  very  accurately  at  the  expense  of  CPU  resources,  or  it 
may  be  approximated  to  a  lower  degree  of  accuracy  with  little  CPU  resources  expended. 

At  this  point,  a  simple,  efficient  approach  for  finding  the  distance  from  the  cell  center 
to  the  closest  wall  was  chosen.  The  algorithm  begins  by  gathering  all  of  the  faces  on  the 
boundary  patches  which  compose  our  solid  surface.  From  this  list  of  boundary  faces,  a 
simple  search  is  done  finding  the  minimum  distance  from  each  face  centroid  to  the  given 
cell  centroid.  Once  the  closest  face  centroid  is  determined,  the  distance  to  the  nodes  of 
this  face  are  computed  and  the  minimum  of  the  distances  is  taken  to  be  our  distance  to 
the  wall.  Graphically,  this  is  shown  in  figure  2. 

In  figure  2,  the  face  centroids  are  shown  as  solid  circles,  the  nodes  of  the  closest  face  as 
empty  circle,  the  further  distances  are  represented  as  dashed  lines,  and  the  closest  distance 
is  represented  by  the  solid  line.  In  this  case,  the  face  centroid  labeled  as  F4  is  the  closest 
face,  and  the  node  labeled  Ni  is  the  minimum  distance  from  the  cell  centroid.  Therefore, 
in  this  simple  algorithm,  this  distance  would  be  used  as  the  distance  from  the  cell  centroid 
to  the  wall  (d). 

Time  Integration 

The  time-dependent  formulation  of  the  governing  equations  has  been  chosen  and  is 
valid  for  time-dependent  and  steady  state  problems.  The  very  small  time  scales  associated 
with  finite-rate  chemistry  can  often  make  completely  explicit  time  integration  strategies 
for  flows  with  strong  chemical  reactions  impractical.  Linearization  of  the  chemical  reaction 
source  terms  helps  to  eliminate  the  stability  constraints  associated  with  the  chemical  reac¬ 
tions  and  allows  the  use  of  practical  time  steps  in  performing  both  steady  state  and  time 
dependent  calculations.  Once  the  linearizations  of  the  various  components  of  the  govern¬ 
ing  equations  have  been  obtained,  they  may  be  coupled  with  any  of  the  current  implicit 
time  integration  techniques  for  unstructured  discretizations.  Thus  any  advances  in  time 
integration  techniques  for  perfect  gas  unstructured  discretions  may  be  quickly  applied  to 
the  current  algorithm. 

One  implicit  time  integration  scheme  which  has  received  a  lot  of  recent  interest  is 
the  Generalized  Minimum  Residual  (GMRES)  algorithm  originally  proposed  by  Saad  and 
Schultz[31].  The  GMRES  algorithm  has  been  successfully  adapted  to  both  structured 
flow  solvers[32]  and  unstructured  flow  solvers [33].  The  GMRES  implicit  time  integration 
technique  is  one  of  several  chosen  for  implementation  in  the  unstructured  flow  solver.  The 
others  include  explicit  m-stage  Runge-Kutta  with  implicit  residual  smoothing  and  Jacobi 
iteration  as  an  inner  iterative  technique  for  the  Euler  implicit  method  [30] . 

Results 

Test  cases  were  sought  which  would  adequately  test  the  individual  components  of  the 
unstructured  fluid  dynamics  code.  Here  we  present  several  of  these  test  cases. 

The  first  of  these  solutions  is  chemically  reacting  air  flowing  through  a  supersonic 
nozzle.  The  air  is  assumed  to  be  in  chemical  and  thermodynamic  equilibrium. 

The  second  case  corresponds  to  the  chemically  reacting  flow  field  about  the  Ramll-C 
reentry  probe.  For  this  problem  the  air  is  assumed  to  be  in  chemical  and  thermodynamic 
non-equilibrium . 
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The  third  solution  presented  is  the  aeroassist  flight  experiment  (AFE)  at  Mach  10  and 
Mach  31.7.  Solutions  for  the  Mach  10  test  case  are  compared  to  experimental  data.  The 
Mach  31.7  solutions  are  compared  to  other  numerical  studies  of  the  AFE. 

The  fourth  solution  presented  is  an  analytic  forebody.  The  analytic  forebody  is  a  three 
dimensional,  inviscid,  supersonic  flow  over  the  cockpit  of  a  high  performance  aircraft.  The 
analytic  forebody  has  been  studied  both  numerically  and  experimentally  [34]  This  problem 
was  chosen  to  test  the  three  dimensional  k-exact  reconstruction  algorithm. 

The  last  set  of  solutions  presented  are  designed  to  test  the  viscous  fluxes  and  the 
turbulence  model.  The  first  solution  in  this  set  is  the  subsonic,  laminar  flow  over  a  two 
dimensional  flat  plate.  The  numerical  solution  is  compared  against  the  analytical  solution 
of  Blasius.  The  second  solution  is  the  supersonic,  turbulent  flow  over  the  same  fiat  plate. 
For  this  case,  solutions  are  compared  against  the  experimental  data  curve  fit  of  Clauser 

[35]. 

Axi-symmetric  Nozzle  ’With  ICang  &  Dunn  Air  Chemistry 

The  first  solution  presented  is  chemically  reacting  air  flowing  through  a  diverging 
supersonic  nozzle.  The  cross  section  of  the  nozzle  is  circular  with  the  cross  sectional  radius 
described  by  the  function 

r  =  0.125  L[l  +  sin(7rA:/2L)] ,  L  =  2m 
where  X  is  the  streamwise  location. 

Using  the  Kang  &  Dunn  air  model  [36],  solutions  were  obtained  for  inflow  conditions 
corresponding  to  air  in  thermodynamic  and  chemical  equilibrium  at  a  Mach  number  of 
1.39  with  a  temperature  of  9000  K  and  a  mixture  density  of  3.856  x  10  ^Kg/m^. 

Solutions  were  obtained  on  three  different  grid  types;  a  structured  hexahedral  grid, 
a  fully  unstructured  hexahedral  grid,  and  an  unstructured  prism  grid.  In  addition,  the 
solution  was  obtained  on  5  different  levels  of  grid  refinement  for  each  mesh  type,  with 
level  1  representing  the  coarsest  mesh  and  level  5  representing  the  finest  mesh.  Previous 
solutions  for  the  unstructured  meshes  were  interpolated  to  the  next  mesh  level  using  a 
simple  quad  tree  algorithm. 

Figures  3(a)  through  3(d)  show  the  level  1  through  level  4  prism  meshes.  Figures 
4(a)  through  4(d)  show  the  level  1  through  level  4  hexahedral  meshes.  Provided  with  the 
figures  is  the  nrunber  of  cells,  faces,  and  nodes  for  each  mesh. 

Figures  5(a)  and  5(b)  show  the  centerline  N2  mass  fraction  profiles  for  the  level  4  and 
level  5  meshes.  From  these  two  figures,  we  see  that  grid  refinement  has  little  effect  on  the 
mass  fraction  profiles  for  this  problem.  In  addition,  the  solutions  on  the  different  mesh 
types  agree.  Also  evident  from  these  figures  is  that  assuming  equilibrium  mass  fractions 
in  the  inflow  has  little  effect  on  the  mass  fraction  profiles. 

Figures  6(a)  and  6(b)  show  the  temperature  profile  along  the  centerline.  From  these 
figures  we  see  a  large  dependency  between  the  different  meshes.  The  prism  solutions  differ 
rather  dramatically  from  the  two  hexahedral  solutions  towards  the  center  of  the  nozzle.  In 
addition,  the  solutions  on  the  different  mesh  levels  differ  considerably.  Figure  6(c)  shows 
the  centerline  temperature  profile  for  all  of  the  prism  meshs,  and  figure  6(d)  shows  the 
profiles  for  all  of  the  unstructured  hexahedral  meshes.  From  figures  6(c)  and  6(d)  we 
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can  clearly  see  the  grid  dependency  of  the  centerline  temperature  profile.  It  is  interesting 
to  note  that  the  prism  solution  is  approaching  a  grid  converged  solution  faster  than  the 
hexahedral  mesh  despite  having  less  cells  on  meshes  of  the  same  level. 

Figures  7(a)  and  7(b)  show  the  N2  mass  fraction  profiles  along  the  surface  of  the  nozzle 
for  level  4  and  level  5  meshes,  respectively.  Again,  we  see  little  effect  of  grid  convergence 
or  the  assumption  of  equilibrium  mass  fractions. 

Figmes  8(a)  and  8(b)  show  the  temperature  profiles  on  the  surface  of  the  nozzle 
for  level  4  and  level  5  meshes,  respectively.  Unlike  the  centerline  profiles  which  showed 
a  strong  grid  dependency  towards  the  center  of  the  nozzle,  the  surface  profiles  show  a 
fairly  strong  grid  dependency  at  the  inflow.  ^Ve  also  can  clearly  see  the  effect  of  assuming 
equilibrium  conditions  at  the  inflow.  As  the  grid  is  refined,  the  inflow  is  tending  away 
from  the  equilibrium  mass  fractions,  as  expected.  Again,  we  see  that  the  prism  meshes  are 
capturing  the  trends  earlier  than  the  hexahedral  solutions. 

Ram  II-c 

During  the  late  sixties  a  series  of  probes  were  flown  in  a  velocity  regime  in  which 
ionization  would  play  an  important  role.  The  second  series  of  tests,  known  as  the  Ram  II-c 
experiments,  have  become  a  popular  test  case  for  modern  CFD  codes.  During  the  later 
series  of  flights,  water  and  electrophilic  liquid  [36]  were  injected  into  the  plasma  layer  to 
reduce  the  free  electron  density  level.  However,  only  data  for  the  flights  without  injection 
are  available. 

The  geometry  of  the  probes  consisted  of  a  sphere-cone  configuration  with  a  15.24 
cm  nose  radius  and  a  9  degree  cone  half  angle.  The  probes  were  instrumented  to  measure 
electron  density  during  the  flight.  Eight  wire  probes  were  located  on  the  leading  edge  of  the 
rake  with  the  last  probe  extending  approximately  7  cm  into  the  plasma  layer.  Collectors 
were  placed  with  their  longitudinal  axis  at  46  degrees  to  the  flow  direction  and  biased  with 
a  constant  negative  voltage  to  collect  ions. 

The  case  considered  in  this  report  corresponds  to  a  flight  altitude  of  61Km  with  free 
stream  conditions;  Too  —  254A’^,  {7oo  —  7650m/s,  and  a  Reynolds  number  based  on  the 
nose  radius  Ren  =  19, 500.  These  conditions  correspond  to  a  flight  Mach  number  of  23.9. 
The  chemistry  model  used  to  obtain  the  solutions  is  the  air  chemistry  model  described  by 
Park  [37].  The  Park  air  chemistry  model  contains  7  species  (A2, 02,  JVO,  ATO"*",  AT, O, e”) 
and  18  reactions. 

Figures  9(a)  and  9(b)  show  the  level  2  and  level  3  prism  meshes  used  to  obtain  the 
unstructured  solution.  The  computed  Mach  contours  for  the  two  meshes  are  shown  in 
figures  10(a)  and  10(b). 

Figure  11  shows  the  electron  number  density  on  the  surface  of  the  probe  for  the 
unstructured  prisms,  structured  hexahedrals  and  experimental  data.  On  these  meshes, 
both  the  structured  hexediedral  and  the  unstructured  prism  solutions  overpredicted  the 
electron  number  density  on  the  aft  portion  of  the  body.  It  is  worth  noting  that  if  the 
flowfield  was  initialized  to  the  freestream  Mach  number  that  convergence  problems  for 
both  the  structured  and  unstructured  codes  were  experienced.  In  this  case,  the  shock  wave 
forms  on  the  body  and  moves  outward  into  the  flow  field  resulting  in  a  severe  transient 
which  limited  the  time  step.  By  initializing  the  flow  to  a  low  Mach  number  and  enforcing 
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freestream  conditions  on  the  outer  boundary,  the  shock  moves  inward  towards  the  body, 
enabling  a  much  higher  Courant  number. 

Aeroassist  Flight  Experiment 

A  family  of  vehicles,  known  as  aeroassist  vehicles,  are  currently  under  design.  These 
vehicles  operate  in  the  upper  reaches  of  the  atmosphere  and  at  velocities  higher  than  those 
encountered  by  typical  reentry  vehicles. 

A  class  of  these  vehicles,  known  as  aeroassist  orbital  transfer  vehicles  (AOTV),  are 
being  designed  to  transfer  payloads  from  a  low  earth  orbit  to  a  high  earth  orbit  and  back. 
In  returning,  the  vehicles  velocity  must  be  greatly  reduced  to  the  earth  orbital  velocity 
and  achieve  capture  in  low  earth  orbit.  This  reduction  in  velocity  may  be  achieved  by  one 
of  two  methods,  the  use  of  retro  rockets  or  by  allowing  aerodynamic  drag  forces  to  act  on 
the  vehicle.  Studies  indicate  that  fuel  loads  are  lowered  and  therefore  the  payload  may 
be  increased  using  the  second  method.  Because  of  the  high  altitude  and  high  velocity, 
ground  based  facilities  are  not  capable  of  reproducing  these  conditions.  In  order  to  better 
understand  the  physics  of  these  problems,  the  AFE  has  been  proposed  by  the  National 
Aeronautics  and  Space  Administration.  The  vehicle  will  be  transported  by  the  space 
shuttle  and  placed  in  a  low  earth  orbit.  The  AFE  will  then  be  propelled  into  the  atmosphere 
to  simulate  the  velocity  and  trajectory  of  a  return  mission  from  geosynchronous  orbit.  The 
data  obtained  will  then  be  used  to  validate  current  computational  fluid  dynamics  codes 
which  will  be  used  in  future  AOTV  designs. 

The  AFE  will  travel  in  the  Earths  upper  atmosphere  at  velocities  ranging  from  7  to 
10  km/s.  At  these  conditions,  chemical  and  thermal  non-equilibrium  effects  are  significant. 
The  AFE  project  was  designed  to  obtain  critical  flight  data  that  will  be  used  to  validate 
hypersonic  computational  fluid  dynamic  codes.  Two  sets  of  solutions  for  the  AFE  were 
obtained.  The  first  set  of  solutions  corresponds  to  a  free  stream  Mach  number  of  10.  The 
Mach  10  case  is  significant  because  good  experimental  data  exists  for  this  problem.  The 
second  set  of  solutions  is  for  a  free  stream  Mach  number  of  31.7.  In  this  case  only  code  on 
code  validations  are  available  at  this  time. 

Solutions  to  the  AFE  proved  to  be  very  stiff  and  difiicult  to  obtain.  The  solutions  pre¬ 
sented  here  were  obtained  with  the  GM-Res  time  integration  scheme.  To  further  increase 
the  efficiency,  mesh  sequencing  was  used  on  a  series  of  meshes.  A  total  of  5  meshes  were 
generated  in  the  process.  Several  views  of  the  flnest  AFE  mesh  are  shown  in  figure  12.  The 
body  of  the  AFE  has  been  mirrored  about  the  plane  of  symmetry  to  provide  a  full  view 
of  the  proposed  configuration.  The  mesh  displayed  in  figure  12  contains  approximately 
106,000  cells  and  202,000  faces. 

The  solution  for  the  Mach  10  test  case  is  presented  in  flgure  13.  This  figure  shows  a 
comparison  of  the  computed  and  experimental  data  for  the  pressure  coefficient  over  the 
surface  of  the  AFE  in  the  plane  of  symmetry.  Two  sets  of  computed  results  are  presented 
in  this  figure;  the  first  set  corresponds  to  a  medium  mesh  of  approximately  50,000  cells, 
and  the  second  set  corresponds  to  the  fine  mesh  shown  in  figure  12.  As  can  be  seen  from 
the  figure,  both  sets  of  computed  results  compare  favorably  to  the  experimental  data  of 
J.R.  Micol  [38]. 

Results  for  the  Mach  31.7  test  case  are  presented  in  figure  14(a)  through  14(d).  The 
freestream  conditions  for  this  test  case  correspond  to  point  1  of  the  planned  trajectory. 
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This  corresponds  to  a  free  stream  velocity  of  9.818  Km/s,  a  mixture  density  of  1.542  x 
10”®  kg/m®,  and  a  temperature  of  195.3  K.  The  flow  was  assumed  to  be  in  chemical  and 
thermal  non-equilibrium.  The  freestream  mass  fractions  were  obtained  using  a  standard 
air  composition  for  N2  and  O2.  All  other  species  were  absent  from  the  freestream.  The 
solutions  were  obtained  using  the  7  species,  18  reaction  Park  air  chemistry  model  [37] 
discussed  earlier. 

In  flgures  14(a)  and  14(b),  the  electron  number  density  and  temperature  are  shown. 
The  electron  number  density  plot  compares  favorably  to  the  solutions  obtained  by  Greendyke, 
Gnoffo,  and  Lawrence  [39].  In  flgures  14(c)  and  14(d),  contours  of  the  O2  and  O  species 
densities  are  shown.  At  these  temperatures  we  would  expect  to  see  complete  dissociation 
of  O2.  This  is  clearly  evident  from  this  figure. 

It  should  be  pointed  out  that  the  assumptions  in  the  thermo- chemical  modeling  im¬ 
plemented  in  the  unstructured  code  are  not  entirely  valid  for  this  flow  field.  In  addition, 
the  grid  for  the  Mach  31.7  case  needs  to  be  refined  near  the  body  of  the  vehicle.  For  the 
simulation  of  the  high  speed  test  case  the  mesh  is  far  too  coarse  behind  the  shock  to  prop¬ 
erly  capture  all  of  the  physics  of  this  complex  flow  field.  However,  despite  the  assumptions 
and  poor  quality  of  the  mesh,  the  results  are  adequate.  The  effect  of  the  assumptions  is 
also  seen  for  the  Ram-IIc  simulations  shown  in  Applebaum  [30]  [40] 

Analytic  Forebody 

The  analytic  forebody  problem  has  been  studied  experimentally  as  well  as  numerically 
[34].  The  problem  represents  the  inviscid  supersonic  flow  over  the  cockpit  region  of  a  high 
performance  aircraft.  The  geometry  consist  of  a  cross  section  which  develops  smoothly 
from  circular  at  the  nose  through  a  lobed  analytic  curve  and  back  to  circular  at  the  aft 
section  of  the  body.  This  case  was  chosen  to  test  the  k-exact  reconstruction  method  in 
three  dimensions.  The  solution  was  obtained  using  mesh  sequencing  on  two  meshes.  The 
finer  of  the  meshes  is  shown  in  figure  15.  This  mesh  corresponds  to  a  structured  31x31x41 
mesh  which  has  been  subdivided  into  tetrahedral  cells.  Three  views  of  the  forebody  are 
shown,  clockwise  from  the  top,  they  represent  the  surface  of  the  forebody,  exit  plane,  and 
symmetry  plane.  The  mesh  in  figure  15  contains  approximately  230,000  cells,  440,000 
faces,  and  40,000  nodes. 

The  freestream  conditions  for  this  problem  are  given  by:  Moo  =  1-7,  Poo  —  1.1774 
kglm^,  Poo  =  101325  N/m^,  and  the  angle  of  attack  a  =  0.  The  solution  was  obtained 
using  the  k-exact  reconstruction  method  for  the  inviscid  fluxes.  The  inviscid  fluxes  were 
computed  using  the  Roe  flux  difference  splitting  algorithm  discussed  earlier.  The  Runge- 
Kutta  explicit  time  integration  scheme  was  used  to  advance  the  solution  in  time.  The 
solution  on  the  coarse  mesh  required  11306  iterations  to  reduce  the  residual  6  orders.  On 
a  Cray  YMP  this  required  17  minutes  of  CPU  time.  The  finer  solution  was  converged  an 
additional  2  orders  in  8306  iterations  which  required  1  hour  and  46  minutes  of  Cray  YMP 
time. 

The  solution  to  the  forebody  is  shown  in  figures  16  and  17.  In  figure  16,  the  pressure 
contours  are  shown  on  the  top  and  bottom  planes  of  symmetry.  In  figure  17,  the  surface 
pressure  distribution  in  the  plane  of  symmetry  is  compared  to  experimental  data  [34] .  As 
can  be  clearly  seen,  the  results  compare  favorably  to  the  experimental  data.  Note  that  in 
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the  aft  section  of  the  forebody  the  predicted  values  are  slightly  lower  than  the  experimental 
data.  Papay  [41]  showed  that  for  an  inviscid  solution  to  the  analytic  forebody  this  is 
typical.  Papay  also  showed  that  solutions  obtained  with  with  an  appropriate  turbulence 
model  provide  better  agreement  in  this  region. 

2-D  Laminar  Flat  Plate 

The  two  dimensional  flow  of  a  perfect  gas  over  a  one  meter  long  flat  plate  at  a 
freestream  Mach  number  of  0.3  is  presented  in  figure  18.  The  freestream  temperature 
and  density  are  500  K  and  0.03973  kg  jm^  respectively.  The  Reynolds  number  based  on 
the  length  of  the  plate  is  200,000  which  is  well  within  the  valid  range  for  laminar  flow.  The 
problem  is  aligned  so  that  the  plate  runs  along  the  x-axis  and  the  y-axis  is  normal  to  the 
plate.  To  obtain  the  solution,  a  quadrilateral  21x41  mesh  was  converted  to  a  triangular 
mesh  by  subdivided  each  quadrilateral  into  two  triangular  cells. 

The  boundary  condition  at  the  inflow  fixes  the  total  temperature  and  total  pressure. 
The  pressure  is  extrapolated  from  the  interior,  and  the  density  and  velocity  are  computed 
based  on  the  extrapolated  pressure  and  the  specified  total  temperature  and  total  pressure. 
The  boundary  condition  at  the  outflow  sets  the  pressure  to  the  specified  back  pressure. 
All  other  flow  quantities  are  extrapolated.  The  surface  of  the  flat  plate  is  modeled  using 
a  no-slip  adiabatic  boundary  condition.  This  boundary  condition  extrapolates  the  species 
densities,  pressure,  and  temperature  from  the  first  interior  cell.  The  velocity  components 
u,v,w  are  set  to  zero.  For  the  top  boundary,  all  values  are  extrapolated  from  the  interior. 

Figure  18  shows  a  comparison  of  the  similarity  profiles  computed  at  the  exit  plane 
to  those  of  Blasius.  Two  solutions  are  presented  in  this  figure.  One  of  the  solutions  was 
computed  using  the  k-exact  reconstruction  method  for  the  viscous  flux.  The  other  solution 
was  obtained  using  the  thin  layer  approximation.  The  solutions  were  obtained  using  the 
GMRES  implicit  time  integration,  along  with  a  second  order  k-exact  reconstruction  for 
the  inviscid  fluxes. 

Both  solutions  provide  excellent  results.  The  numerical  boundary  layer  thickness  is 
approximately  0.012  meters  which  compares  favorably  to  the  Blasius  value  of  0.0116. 

2-D  Turbulent  Flat  Plate 

The  solution  presented  in  this  section  represents  the  two  dimensional  turbulent,  su¬ 
personic  flow  over  a  one  meter  long  flat  plate.  The  conditions  in  the  freestream  correspond 
to  Moo  =  2.244,  Too  =  170K,  and  poo  =  0A538kg/m^.  This  corresponds  to  a  Reynolds 
number  of  2x10'’’.  The  solutions  were  obtained  using  the  one  equation  Spalart  and  All- 
maras  model  discussed  earlier.  The  turbulence  model  working  variable  was  initialized  to 
ten  times  the  freestream  laminar  viscosity.  The  k-exact  reconstruction  algorithm  was  used 
to  determine  the  viscous  and  inviscid  fluxes.  The  wall  boundary  condition  extrapolates 
the  density,  pressure,  and  temperature  and  sets  the  velocity  components  to  zero.  This 
corresponds  to  an  adiabatic  wall  boundary.  To  obtain  the  solution,  a  41x81  quadrilateral 
mesh  was  converted  to  a  triangular  mesh.  This  was  accomplished  by  subdividing  each 
quadrilateral  cell  into  two  triangular  cells.  The  solution  was  advanced  in  time  using  the 
Runge-Kutta  explicit  time  integration  scheme.  The  solution  required  59853  iterations  to 
converge  the  solution  6  orders.  On  a  Cray  YMP  this  required  approximately  23  minutes. 
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Figure  19  compares  the  unstructured  solution  at  the  exit  plane  to  the  analytic  law 
of  the  wall  plot.  The  analytical  solution  was  obtained  using  the  coefficients  developed  by 
Clausei-  [35].  Using  these  coefficients,  the  equation  for  the  log  region  is  given  by 

u+ =  5.61og(j/+)  +  4.9 


In  this  figure,  the  frictional  velocity  u*  was  taken  to  be  constant  throughout  the  boundary 
layer.  The  value  of  the  frictional  velocity  was  obtained  from 


V 


The  wall  shear  stress  was  evaluated  using  the  vorticity  at  the  first  cell  off  the  wall. 


Conclusions 

A  method  for  the  solution  to  the  Navier  Stokes  equations  with  a  generalized  thermo¬ 
chemical  model  and  one  equation  turbulence  model  on  unstructured  meshes  was  presented. 
In  addition  a  simple  method  for  the  k-exact  reconstruction  was  reviewed.  Several  test  cases 
were  presented  to  test  the  efficiency  of  the  viscous  and  thermo-chemical  modeling.  From 
the  test  cases,  it  was  seen  that  the  current  methodology  produces  accurate  results  for  the 
solutions  presented. 


Future  Work 

The  aim  of  our  generalized  solution  development  program  is  to  be  able  to  accurately 
and  efficiently  solve  the  Navier  Stokes  equations  on  unstructured  meshes  with  both  equilib¬ 
rium  and  non-equilibrium  chemistry  and  thermodynamics.  Building  on  this  goal  it  is  our 
hope  to  implement  more  turbulence  models  into  the  unstructured  flow  solver.  Extensive 
validation  of  the  flow  solver  is  required. 
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Figure  3(a):  Air  nozzle  level  01  prism  mesh 
250  cells,  271  faces,  106  nodes 


Figure  3(b):  Air  nozzle  level  02  prism  mesh 
958  cells,  1077  faces,  362  nodes 


Figure  3(c):  Air  nozzle  level  03  prism  mesh 
3742  cells,  4287  faces,  1328  nodes 


Figure  3(d):  Air  nozzle  level  04  prism  mesh 
14798  cells,  17109  faces,  5090  nodes 


Figiire  4(a):  Air  nozzle  level  01  hexahedral  mesh 
350  cells,  425  faces,  252  nodes 


Figure  4(b):  Air  nozzle  level  02  hexahedral  mesh 
1300  cells,  1650  faces,  902  nodes 


Figure  4(c):  Air  nozzle  level  03  hexahedral  mesh 
5000  cells,  6500  faces,  3402  nodes 


Figure  4(d):  Air  nozzle  level  04  hexahedral  mesh 
19600  cells,  25800  faces,  13202  nodes 


N,  mass  fraction  N,  mass  fraction 


Figure  5(a):  Air  nozzle  centerline  N2  mass  fraction 
profile  for  level  04  grids 


Figure  5(b):  Air  nozzle  centerline  N2  mass  fraction 
profile  for  level  05  grids 
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Figure  9(b):  Ram  II-C  level  03  prism  mesh 
21686  cells,  25145  faces,  7386  nodes 
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Figure  10(b):  Ram  II-C  Mach  contours  on  level  03  mesh 
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Figure  11:  Riim  II-C  surface  electron  number  density 

at  61Km 
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Figure  12:  Several  views  of  the  Aeroassist  Flight  Experiment  mesh. 


Fig  13:  Pressure  distribution  on  surface 
of  AFE  for  Mach  10  test  case 
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Fig  17:  Analytic  Forebody  surface  pressure  distribution 
in  symmetry  plane 
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Fig  18:  Similarity  profiles  for  subsonic 
laminar  fiat  plate 


Fig  19:  Wall  law  plot  for  the  supersonic 
flat  plate 


